home *** CD-ROM | disk | FTP | other *** search
/ PsL Monthly 1993 December / PSL Monthly Shareware CD-ROM (December 1993).iso / prgmming / dos / c / hzip.com / PKMG.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1991-01-08  |  7.0 KB  |  250 lines
  1. /////////////////////////////////////////////////////////////////////
  2. // pkmg.cpp: Length restricted huffman code generator routines.
  3. // Copyright (c) 1991 Azarona Software
  4. // All rights reserved. 
  5. //
  6. // Implements the package-merge algorithm as described in the
  7. // paper "A Fast Algorithm for Optimal Length-Limited Huffman
  8. // Codes", Lawrence L. Larmore and Daniel S. Hirschberg, JACM, 
  9. // Vol. 37, No.3, July 1990, pp. 464-473
  10. //
  11. // This code is my interpretation of the algorithm, (it's presented
  12. // rather abstractly in the paper). I hope I got it right.
  13. //
  14. // The straight-forward implementation of the algorithm consumes
  15. // lots of memory. This implemetation has been optimized as much
  16. // as I know how, reducing the memory usage from over 200k, to 
  17. // around 64k. The authors Larmore and Hirschberg describe a way
  18. // to use much less space, but they don't give the full algorithm
  19. // and I couldn't make heads or tails out of it.
  20. //
  21. // -- Bryan Flamig
  22. /////////////////////////////////////////////////////////////////////
  23.  
  24. #include <stdio.h>
  25. #include <stdlib.h>
  26. #include <alloc.h>
  27. #include "pkmg.h"
  28.  
  29. void record_deepest(pnode *p, char soln[])
  30. // soln[sym] will eventually contain the leaf depth for
  31. // symbol sym. So, for all symbols in the package p, we
  32. // we record the deepest level found for those symbols.
  33. {
  34.   while (p) { // Go thru complete package
  35.     if (p->lvl != -1) { 
  36.       // We have a leaf pnode, might replace level in solution set
  37.       if (p->lvl > soln[p->sym]) soln[p->sym] = p->lvl;
  38.     }
  39.     p = p->n; // Next pnode of package
  40.   }
  41. }
  42.  
  43. int next_diadic(unsigned n)
  44. // Returns the power of two exponent for the first bit
  45. // in n having a 1 in it, scanning right to left
  46. {
  47.   int p = 0;
  48.   if (n == 0) { printf("Error: infinite diadic\n"); exit(1); }
  49.   while(1) {
  50.     if (n & 1) return p;
  51.     n >>= 1;
  52.     p++;
  53.   }
  54. }
  55.  
  56. int package_merge(long wts[], unsigned char smap[], 
  57.                   int ns, int nl, char soln[])
  58. // The weigts wts[] must come in sorted in decreasing order. If 
  59. // you don't do this, all bets are off. (The program might even
  60. // crash!) smap[] = the indirection array, where smap[i] is the
  61. // symbol value for the ith element. ns = total number of symbols
  62. // to work with, nl = max number of levels you're going to allow, 
  63. // soln[] is the list of leaf depths for the symbols.
  64. // NOTE: keep ns * nl <= 4096.
  65. {
  66.   // We alternate between these two nodesets, each of which
  67.   // represents the candidates of a particular level of the
  68.   // optimum nodeset. The packages on the nodesets are always
  69.   // in order, smallest weight package first.
  70.  
  71.   nodeset c1(ns); 
  72.   nodeset c2(ns); 
  73.  
  74.   // The current nodeset we're using
  75.  
  76.   nodeset *candidates;
  77.  
  78.   // X is the value that the width
  79.   // of the optimum nodeset must be
  80.  
  81.   unsigned X = ns - 1;
  82.  
  83.   // Okay, now we have current level, current width (representing
  84.   // only the power of two exponent of the width), the current
  85.   // minimum width, and the success return value for the algorithm
  86.  
  87.   int level, width, min_width, rv;
  88.  
  89.   // Allocate room for nodesets
  90.  
  91. #ifdef DEBUG
  92.   printf("bytesleft = %lu\n", coreleft());
  93. #endif
  94.  
  95.   pnode::pal = new pnode_allocator(ns * nl + 2*ns);
  96.  
  97. #ifdef DEBUG
  98.   printf("nodesleft = %d\n", pnode::pal->room());
  99.   printf("bytesleft = %lu\n", coreleft());
  100. #endif
  101.  
  102.   // Start with the deepest level
  103.  
  104.   level = nl;
  105.   width = -level; // Just using the power of 2 exponent
  106.  
  107.   // Get minimum width from the next diadic value of X, 
  108.   // smallest comes first
  109.  
  110.   min_width = next_diadic(X); // Just using power of 2 exponent
  111.  
  112.   // Okay, initialize the nodeset for this level
  113.   // (assumes c2 = empty nodeset)
  114.  
  115.   c1.setup_and_merge(c2, wts, smap, ns, level);
  116.   candidates = &c1;
  117.  
  118.   while (X > 0) {
  119.  
  120. #ifdef DEBUG
  121.     printf("nodesleft = %d\n", pnode::pal->room());
  122. #endif
  123.  
  124.     if (candidates->isempty() || width > min_width) {
  125.        rv = 0;
  126.        goto get_out;
  127.     }
  128.  
  129.     if (width == min_width) {
  130.        // Delete smallest weight package from candidate nodeset, add to
  131.        // to the solution nodeset
  132.        pnode *p = candidates->remove_front();
  133.        record_deepest(p, soln);
  134.        pnode::rmv_package(p);  // Delete the whole package
  135.        X &= ~(1 << min_width); // Assumes min_width non-negative
  136.        if (X == 0) {           // We have solution!
  137.           rv = 1;
  138.           goto get_out;
  139.        }
  140.        min_width = next_diadic(X);
  141.     }
  142.  
  143.     // Bundle up the candidates into packages, and then
  144.     // merge them with the next highest level candidates
  145.  
  146.     candidates->package(); 
  147.     if (level > 1) { // If at level 0, don't need to merge anymore
  148.        if (candidates == &c1) { // Alternate between nodesets
  149.           c2.setup_and_merge(*candidates, wts, smap, ns, level-1);
  150.           candidates = &c2;
  151.        }
  152.        else {
  153.           c1.setup_and_merge(*candidates, wts, smap, ns, level-1);
  154.           candidates = &c1;
  155.        }
  156.     }
  157.     width++;
  158.     if (level > 0) level--;
  159.   }
  160. get_out:
  161.   delete pnode::pal;
  162.   return rv;
  163. }
  164.  
  165. void make_codes(char len[], unsigned code[], unsigned char smap[],
  166.                 int ns, int nl)
  167. // Given the leaf depths len[], (which also happen to be the 
  168. // code lengths), compute the code for each symbol
  169. {
  170.   int    i;
  171.  
  172.   unsigned *start = new unsigned[nl+2];
  173.   int *len_cnt = new int[nl+1];
  174.  
  175.   // Compute the number of codes at each level
  176.  
  177.   for (i=1; i<nl+1; i++) len_cnt[i] = 0;
  178.   for (i = 0; i<ns; i++) len_cnt[len[i]]++; 
  179.  
  180.   // Compute the starting code for each level
  181.  
  182.   start[1] = 0;
  183.   for (i = 1; i <= 16; i++)
  184.       start[i + 1] = (start[i] + len_cnt[i]) << 1;
  185.  
  186.   // Compute the code for each symbol, mapping back to
  187.   // proper indices while we're at it
  188.  
  189.   for (i = 0; i < ns; i++) 
  190.       code[smap[i]] = (start[len[i]]++) << (16-len[i]);
  191.  
  192.   delete start;
  193.   delete len_cnt;
  194. }
  195.  
  196. void qs(long wts[], unsigned char smap[], int n)
  197. // This specialized quick sort is 60 times faster than
  198. // Turbo C++'s generic qsort() for n = 256
  199. {
  200.   int i, j, l, r;
  201.   long tw, vw;
  202.   unsigned char ts;
  203.   int *top, *stk;
  204.  
  205.   stk = new int[n]; // Actually, only approx. lgN needed 
  206.  
  207.   l = 0; r = n-1; top = stk;
  208.   *top++ = l;
  209.   *top++ = r;
  210.  
  211.   do {
  212.     if (r > l) {
  213.  
  214.        // begin: i = partition(l, r);
  215.        vw = wts[r]; i = l-1; j = r;
  216.        while(1) {
  217.          do { i++; } while(wts[i] > vw);
  218.          do { j--; } while(wts[j] < vw);
  219.          tw = wts[i]; ts = smap[i];
  220.          if (j >= i) {
  221.             wts[i]  = wts[j];  wts[j]  = tw;
  222.             smap[i] = smap[j]; smap[j] = ts;
  223.          }
  224.          else break;
  225.        }
  226.        wts[i]  = wts[r];  wts[r] = tw;
  227.        smap[i] = smap[r]; smap[r] = ts;
  228.        // end: i = partition(l, r)
  229.  
  230.        if ((j)>(r-i)) {
  231.           *top++ = l;
  232.           *top++ = j;
  233.           l = i+1;
  234.        }
  235.        else {
  236.           *top++ = i+1;
  237.           *top++ = r;
  238.           r = j;
  239.        }
  240.     }       
  241.     else {
  242.        r = *(--top);
  243.        l = *(--top);
  244.     }
  245.   } while(top != stk);
  246.  
  247.   delete stk;
  248. } 
  249.  
  250.